# Import the transect data
UF1 <- read.csv("N:/RStor/dmcgrath/Glacier/Cameron Peak Fire/2022 Field Season/Data Analysis/GIS/Distributed Snow/UF1_cluster_slope_aspect.csv") %>%
  dplyr::select(-c(1,2,4:6,9:13,15,17)) %>%
  mutate(Date_Time = as.POSIXct(Date_Time, format = "%m/%d/%Y %H:%M"), date = date(Date_Time), site = "UF1")

BF1 <- read.csv("N:/RStor/dmcgrath/Glacier/Cameron Peak Fire/2022 Field Season/Data Analysis/GIS/Distributed Snow/BF1_cluster_slope_aspect.csv") %>%
  dplyr::select(-c(1,2,4:6,9:13,15,17)) %>%
  mutate(Date_Time = as.POSIXct(Date_Time, format = "%m/%d/%Y %H:%M"), date = date(Date_Time), site = "BF1")

BF2 <- read.csv("N:/RStor/dmcgrath/Glacier/Cameron Peak Fire/2022 Field Season/Data Analysis/GIS/Distributed Snow/BF2_cluster_slope_aspect.csv") %>%
  dplyr::select(-c(1,2,4:6,9:13,15,17)) %>%
  mutate(Date_Time = as.POSIXct(Date_Time, format = "%m/%d/%Y %H:%M"), date = date(Date_Time), site = "BF2")

SNOWEX <- read.csv("N:/RStor/dmcgrath/Glacier/Cameron Peak Fire/2022 Field Season/Data Analysis/GIS/Distributed Snow/SnowEx_cluster_slope_aspect.csv") %>%
  dplyr::select(-c(1,2,4:6,9:13,15,17)) %>%
  mutate(Date_Time = as.POSIXct(Date_Time, format = "%m/%d/%Y %H:%M"), date = date(Date_Time), site = "SNOWEX")

all <- rbind(UF1, BF1, BF2, SNOWEX)
# Import raster data 
bf_sev <- raster("N:/RStor/dmcgrath/Glacier/Cameron Peak Fire/GIS/dnbr_cameronpeak_scaled.tif")
bf_SBS <- raster("N:/RStor/dmcgrath/Glacier/Wyatt/Cameron Peak Fire/GIS/GIS Reference Files/CameronPeak_SBS_final.tif")

#plot(bf_SBS)
xy <- st_as_sf(all, coords = c('Longitude', 'Latitude'), crs = 4326)
SBS <- extract(bf_SBS, xy)
## Warning in .local(x, y, ...): Transforming SpatialPoints to the crs of the
## Raster
combined_all <- cbind(all, SBS) %>%
  mutate(SBS = ifelse(SBS == 0, 1,SBS))
all_daily <- combined_all %>%
  group_by(date, CLUSTER_ID, site) %>%
  summarize(lat = mean(Latitude), long = mean(Longitude), depth = mean(Depth),
            slope = mean(slope), aspect = mean(aspect), burn = mfv(SBS)) %>%
  mutate(aspect_dir = ifelse(aspect < 135 & aspect >= 45, "East",
                                    ifelse(aspect < 225 & aspect >= 135, "South",
                                           ifelse(aspect < 315 & aspect >= 225, "West", "North"))),
         burn_class = ifelse(burn >= 2, "Burned", "Unburned"))
## `summarise()` has grouped output by 'date', 'CLUSTER_ID', 'site'. You can
## override using the `.groups` argument.
SD_plots <- ggplot(all_daily, aes(x = aspect_dir, y= depth, color = burn_class)) +
  geom_boxplot(position = "dodge") + 
  facet_wrap(~date) + 
  ggtitle("Transect Snow Depth") + 
  theme(legend.position = "top")

ggplotly(SD_plots)
?geom_boxplot
## starting httpd help server ... done